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Abstract 

In this paper we develop adaptive numerical schemes for certain nonlinear variational prob¬ 
lems. The discretization of the variational problems is done by a suitable frame decomposition 
of the solution, i.e., a complete, stable, and redundant expansion. The discretization yields an 
equivalent nonlinear problem on £2 (Af), the space of frame coefficients. The discrete problem is 
then adaptively solved using approximated nested fixed point and Richardson type iterations. 

We investigate the convergence, stability, and optimal complexity of the scheme. This con¬ 
stitutes a theoretical advantage, for example, with respect to adaptive finite element schemes 
for which convergence and complexity results are usually hard to prove. The use of frames is 
further motivated by their redundancy, which, at least numerically, has been shown to improve 
the conditioning of the corresponding discretization matrices. Also frames are usually easier to 
construct than Riesz bases. Finally, we show how to apply the adaptive scheme we propose for 
finding an approximation to the solution of the PDE governing magnetohydrodynamic (MHD) 
flows, once suitable frames are constructed. 

AMS subject classification: 41A46, 42C14, 42C40, 46E35, 65J15, 65N12, 65N99, 76D03, 76D05, 
76W05 

Key Words: Magnetohydrodynamics, Nonlinear operator equations, Multiscale methods, Overlap¬ 
ping domain decomposition, Adaptive numerical schemes, Frames, Wavelets and multiscale bases. 


1 Introduction 

Adaptive numerical methods have yielded very promising results [2, 3, 4, 6, 8, 9, 12, 22, 47] when 
applied to a large class of operator equations, in particular, PDE and integral equations. In classical 
schemes the adaptivity is realized at the level of the discretization and in the finite element space. 
The finite element space is refined and enriched locally at each iteration step depending on some (a 
posteriori) error estimators [23, 37]. A novel paradigm for adaptive schemes has been recently pro¬ 
posed by Cohen, Dahrnen, and DeVore in [8, 9], where the discretization via wavelet decompositions 
is fixed at the beginning. The adaptivity is indeed realized at the level of the solver of the equivalent 
bi-infinite system of linear equations. The basic idea is to transform the original problem of PDE 
into a discrete (bi-infinite) linear problem on £ 2 (A/"), the space of wavelet coefficients. The discrete 
problem is then solved with the help of approximate iterative schemes. The advantage of the latter 
approach is the fact that its convergence and stability can be proved and its optimal complexity 
can be estimated asymptotically in terms of the number of algebraic operations needed. On the 
contrary, it has been a very hard technical problem to obtain such nice theoretical estimates for 
classical finite element methods, although some important theoretical results has recently appeared 
in [5, 41] for linear elliptic equations. A version of the paradigm in [8, 9] has been recently proposed 
also for nonlinear problems in [10]. It is again based on (wavelet) bases discretizations. 

One drawback of the wavelet approach is the construction of the wavelet system itself especially 
on domains with complicated geometry or manifolds [16, 17, 18]. The wavelet bases constructed so 
far exhibit relatively high condition numbers and limited smoothness. In particular, the patching 
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used to construct global smooth wavelets by domain decomposition techniques appears complicated 
and, in most cases, makes the conditioning even worse. In fact the global smoothness of the basis, 
when implementing adaptive schemes in [8, 9], is a necessary condition for getting compressibility 
(i.e., finitely banded approximations) of (bi-infinite) stiffness matrices, especially for high order op¬ 
erators. This bottleneck has led to generalizations of the Cohen, Dahmen and DeVore approach. 
These generalizations are based on frame discretizations , i.e., stable, redundant, non-orthogonal 
expansions [13, 40], which are much more flexible and simpler to construct even on domains of 
complicated geometry. Frame construction is usually implemented by Overlapping Domain Decom¬ 
positions (ODD) so that patching at the interfaces is no more needed to obtain global smoothness. 
Moreover, the use of frames, due to their intrinsic redundancy, improves the conditioning of the 
corresponding discretization matrices. Certainly, ODD generates regions of the domain where the 
side effect of the redundancy is that functions are no longer uniquely representable by the global 
frame system. At first sight, it may seem that redundancy contradicts the minimality requirement 
on the amount of information being used to approximate the solution. Especially in fluid mechanics, 
accurate simulations already require processing of huge amounts of data. How can one attempt such 
computations if the information is also made redundant? A figurative answer to this question is the 
so-called “dictionary example”: The larger and richer is my dictionary the shorter are the phrases 
I compose. The use of proper terminology avoids long circumlocutions for describing an object. Of 
course, the key point is the capability to choose the right terminology. Back to mathematical terms, 
the combination of adaptivity (i.e., the capability to choose the right terminology) and redundancy 
(i.e., the richness or non-uniqueness of representations) indeed gives rise to fast and accurate approx¬ 
imations [20, 27, 33, 42, 43]. Numerical experiments in [48] show that frames improve conditioning 
without increasing the effective dimension of the problem, i.e. without increasing the number of the 
relevant quantities needed for computations. 

This encourages us to present a generalization of the approach proposed in [10] to (wavelet) frame 
discretizations for some specific nonlinear PDE, in particular those describing certain magnetohy¬ 
drodynamics problems. Magnetohydrodynamics [7, 32, 34, 35, 36, 38, 39, 49] studies macroscopic 
interactions between magnetic field and fluid conductors of electricity. In particular, the following 
physical phenomena are studied: A flow of an electrically conducting fluid across magnetic lines 
causes an electric current in the fluid. The electric current alters the electromagnetic state of the 
system modifying the total magnetic field, which creates the current. The flow of electric current 
across magnetic lines is associated with a body force - Lorentz force - which influences the fluid flow. 
To model the behavior of electrically conducting fluid, the stationary, incompressible Navier-Stokes 
and Maxwell equations, coupled via Ohm’s law and Lorentz force are being used. We refer to [35, 36] 
for a rigorous analysis of a non-adaptive finite-element scheme for the simulation of MHD flows con¬ 
fined to a cubic domain only and arising during electromagnetic purification of molten metals before 
the casting stage. 

We present a way for transforming the nonlinear magnetohydrodynamics problem, possibly defined 
even on more general domains, into an equivalent nonlinear discrete problem on £2 (A/") by using 
suitable frame expansions. We show how the discrete problem can be solved adaptively by means of 
nested fixed point and approximated Richardson iterations. We also discuss convergence, stability, 
and, under certain additional assumptions, computational cost (quasi-optimal complexity) of the 
proposed adaptive procedure. 

The paper is organized as follows: in Section 2 we recall the mathematical model and the 
equations governing MHD flows, specify the boundary conditions, the solution and test function 
spaces. In Section 3, the corresponding weak formulation of the physical problem is derived. Its 
equivalence to a variational nonlinear problem on a suitable subspace of the solution function space 
follows from the standard LBB (Ladyzhenskaya-Babuska Brezzi) theory. The proofs of the results 
listed in this section can be found in [7, 34, 36]. In Section 4, the nonlinear variational problem is 
re-formulated as an equivalent fixed point iteration scheme, where, at each iteration step, a linear 
(non-symmetric) elliptic operator equation is to be solved. Next, we present the way for discretizing 
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the fixed point iteration associated to an abstract nonlinear problem arising in MHD. We translate 
the original nonlinear variational problem into a problem on (.2 (A/ - ), the space of suitable frame 
coefficients. The concept of frames, i.e., stable, redundant, and complete expansions, is recalled in 
Subsection 4.1. In Subsection 4.2 we show how a linear (non-symmetric) elliptic operator equation 
is discretized by means of frame expansions and how Algorithm 1 is used to approximate its solution 
adaptively, up to any prescribed accuracy. Using Algorithm 1 as the main building block, we 
formulate in Section 5 a fully discrete and finite version of the fixed point iteration. We show that 
this discrete version of the fixed point iteration converges to some frame coefficients of the true 
solution of the original problem. In Section 6 we discuss under which conditions on the building 
block procedures in Algorithm 1 the suggested scheme performs quasi-optimally with respect to 
suitable sparseness classes of frame coefficients. In particular, we show how the flexibility and 
redundancy of frames lead to technical difficulties, which do not arise in case of Riesz bases, when 
showing complexity estimates. In Section 7, we show that the relevant solution space for the MHD 
problem is a product of suitable divergence-free vector spaces. We shortly discuss the existence and 
the construction of suitable divergence-free frames for such spaces. This allows for applications of 
the general adaptive scheme we propose. 

Throughout this paper ‘a ~ b' means that both quantities are uniformly bounded by some 
constant multiple of each other. Likewise, ‘a < b' means that there exists a positive constant C such 
that a < Cb. We determine the constants explicitly only if their value is crucial for further analysis. 
The symbol || • ||, when applied to bounded operators, denotes the operator norm from its domain 
space to its image space, these are not always explicitly specified for notational simplicity. 


2 Equations and Boundary Conditions 


In this section, we model the stationary flow of a viscous, incompressible, electrically conducting 
fluid occupying a 3-dimensional bounded region fl C M 3 . For later analysis, we assume that O is a 
Lipschitz domain. Assume also that F, a body force, and E, an externally generated electric field 
are given. Denote by rj the viscosity, p the density, er -1 the electrical resistivity and /./, the magnetic 
permeability of the fluid (some positive constants). 

To describe the interaction of the magnetic field and the electrically conducting fluid in fl math¬ 
ematically we combine the equations of the fluid dynamics and electromagnetic field equations. We 
use Navier-Stokes equations to model the fluid flow and include Lorenz force J x B to express the 
influence of the flow of electric current across magnetic lines on the fluid motion 

—?yAu + p(u • V)u + Vp — J x B = F. (1) 


The flow of the electrically conducting fluid across magnetic lines causes an electric current in the 
fluid. The electric current alters the electromagnetic state of the system modifying the total magnetic 
field, which creates the current in the fluid. This phenomenon is expressed by means of Ohm’s law 

a~ 1 J + \7(/)-u x B = E. (2) 


The total magnetic field B = B 0 + B(J) is decomposed into a sum of a given externally generated 
magnetic field B 0 and the induced magnetic field B(J), which is induced in the fluid by the electric 
current J caused by B. B(J ) is a unique solution (see Lemma 2.2. in [34]) of the quasi-stationary 
form of Maxwell’s equations 


V x p~ 1 B(J) = J and V • B(J) = 0. 
It is also given by the Biot-Savart law 


U(J)(x) = -V x { £ 


4tt Jn l x - y 


J(y U] = -A 


x - y 


4?r J n |x - y| 3 


x J(y)dy, 


(3) 
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for x € R 3 . Note that using (3), we are able to eliminate B from (l)-(2) and solve for J instead. 
Solving for B we would be facing the problem that the magnetic field is defined on all space and 
satisfies different equations inside and outside the region f h The Navier-Stoke’s equations are posed 
inside the region occupied by the fluid, whereas the Maxwell’s equations have to be solved in all of 
space. Thus, the boundary conditions for B on the surface of the region must be specified, which 
is possible only for perfect conductors or perfect insulators. This makes the prescribed boundary 
conditions somewhat artificial. 

We also consider the continuity equations 

V • u = 0, V • J = 0 (4) 

to describe the incompressibility of the fluid and preservation of charge. 

Denote by n the outward unit normal vector field on the boundary of and the stress tensor by 

T(u,p) := — pi + (Vu + (Vu) T ) =: —pi + r/D( u). 

Let the boundary, denoted by T, of the domain Q consist of several relatively open and pairwise 
disjoint components Tj’s, i= 1,... ,4, i.e. T = Id uf 2 uf 3 UF 4 . Different boundary conditions are 
prescribed on each component. The boundary conditions we consider are 
Dirichlet type (prescribed velocity) 

u = gx on T 1; (5) 

Neumann type (prescribed stress) 

Tn = h 2 on T 2 ( 6 ) 


and mixed type for velocity and stress 

u ■ n = g 3 and (Tn -((nT) ■ n)n) =h 3 on T 3 , 

u — (u • n)n = g 4 and (Tn) • n = ft .4 on T 4 . 


(7) 


These boundary conditions are helpful in modeling the free boundary value problems, when dealing 
with the artificially truncated computational domains and the boundary conditions on the artificial 
boundaries. 

Let also the boundary T consist of two relatively open and pairwise disjoint components IF and 
£ 2 , and consider 

Neumann type (prescribed electric current flux through the walls) boundary condition 

J • n = j on £4 ( 8 ) 


and Dirichlet type (prescribed electric potential) boundary condition 


<f> = k on £ 2 . 


(9) 


This helps us to model two different cases: the external magnetic field is given and £4 is not elec¬ 
trically insulating then j ^ 0 ; the magnetic field is generated by the external conductors embedded 
into £ 2 . Note that incorporating the electric potential is useful for various control problems. 


2.1 Function spaces 

Denote by H S (Q) the Sobolev space of square integrable functions v on fi with square integrable 
distributional derivatives D a v up to order s with the norm 

IHI^(0) = (E H^iiLm)’- 

|a| <s 
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For vector-valued functions v = (vi, u 2 , u 3 ) define 


H s (fl) := {v : Vi £ H s (n), i = 1... 3}, 

L 2 (0) := {v : Vi £ L 2 (h),i = 1... 3 

Let H 1 / 2 ( F) denote the fractional order Sobolev space and its dual H 1 ^ 2 ( F)'. 

The minimal regularity assumptions on the given data are 

Fetffll)' EelaiSi), BoetffO), 

g^H 1 / 2 ^), gsGH 1 / 2 (F 3 ), g4 eH 1 / 2 (r 4 ) with g 4 -n = 0 on r 4 , 

h 2 G h 3 £ H 1 / 2 ^)' /j 4 €iL 1/2 ( r 4 )' with h 3 -n = 0 on r 3 . 

Note that (4) imposes the following compatibility conditions on the boundary data 

f gi-n+ f 53 = 0, if F 2 U r 4 = 0, and / j = 0, if S 2 = 0. (10) 

iri J r 3 dsj 

The solution of (l)-(9) is not unique ifr 2 UF 4 = 0 and/or S 2 = 0. This is due to the fact that 
(l)-(9) then only involve the derivatives of p and/or </. To avoid the non-uniqueness of the solution 
we solve for p £ L 2 (Q) (subspace of L 2 (fl) consisting of functions with mean zero) and <fr £ H 1 (fi) 
(subspace of i/ 1 (fi) consisting of functions with mean zero). Therefore, we seek the solution 

u £ H 1 (fl) satisfying (l)-(7), 


J £ L 2 satisfying J • n = j on S 4 
p £ M p := 


L 2 {ft) , ifr 2 ur 4 = 0, 

L 2 ( O) , otherwise 


£ Ma 


H\n) , if e 2 = 0, 

{(j) £ H 1 (VL) : (j) satisfying (9)} , otherwise. 


To derive an equivalent to (l)-(9) variational formulation, we choose the test functions v in the 
space 

Hf(f2) := {v £ H 1 (f2) : v| Fl = 0, v • n|r 3 =0, v - (v • n)n|r 4 = 0}, 

K £ L 2 (fl), q £ M q := M p and 


:= 


H\Q) , if S 2 = 0, 

{q £ 7L 1 (n) : q\s 2 = 0} , otherwise. 


To simplify the notation we introduce the product spaces 

X (U)J) := x L 2 (Q), X (v>Jt) := Hf(fi) x L 2 ( ft), 

:= M p x and M( g> ^) := M q x M$. 


3 Weak Formulation 

In this section we shortly recall how to derive a variational problem equivalent to (1)—(9), (see [7] 
for details). 
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Define a bilinear form ao : X (u j\ xX( u —> R, a trilinear form a 4 :X( u jjxX( u j)XX( u j) —■» R, 

a bilinear form b : X( u j) x (L 2 (fl) x B^Q)) —> R and a linear form x : X( u j) —> R by 


o 0 ((v 1 ,if 1 ),(v 2 ,.If 2 )) := J [ B(v 4 ) : B(v 2 ) + a ” 1 [ if 4 if 2 

1 J o Jo 

+ [ {{K-2 x -Bo) • vi — (Ifi x Bo) • V2), 

Jo 

ai((vi,/fi), (v 2 ,if 2 ), (v 3 ,if 3 )) -=P ((vr • V)v 2 ) • v 3 

Jo 

'3 x £(if 1)) • V 2 - (if2 x B(Jfr)) • v 3 ), 


b((v,K)(q,i/;)) := - f (V • v)g + f K • (W>) 
Jo Jo 


and 

X(v, if) = / F • v + I E ■ K + / h 2 • v + / h 3 • v + / h 4 n • v. 

J ^2 «/ •/r2 *^ 1^3 «/r4 

Let v G Hp(f2) and if G B 2 (^) be arbitrary test functions. Multiplying (1), (2) and (4) by 
corresponding test functions, integrating by parts and using some of the boundary conditions on u 
and J we get the following equivalent variational problem. 

Problem 1: Find (u, J) € X (U)J ) s.f/i. u| Fl = g 4 , (u • n)|r 3 = 93 , (u - (u • n)n)|r 4 = g 4 and 
(p, 0) G such that 

a 0 ((u,J),(v,K)) + a 4 ((u, J), (u, J), (v, if)) + 6((v, if ), (p, 0)) =x(v,if), 

and 

b{(u,J),(q,i>)) = [ jip 
J Si 

/or all (v, if) G X (v , x) and (q,if>) G M {q ^y 


Next, we list some properties (proved in [34, Lemma 2.2] and [7, Corollary 1]) of ao, a 4 and 6 . 

Lemma 3.1. 

a) The forms a 0 , a lr b are bounded on X (u>J ) x X( U)J) , X (u j)XX (u j)XX (u ^ andX( UiJ) xMy )( 
respectively, with 


||ai|| < nrax{p,/u}. 

b) If the intersection o/Hp(f2) and i/ie null space N(D) of the deformation tensor D is {0}, f/ien 
f/ie /orm a 0 is positive definite on X( v x X( v xq *-e. 

oo((v, if), (v, if)) > a 0 ||(v, if)||x (v , JO 

with ao := c(f2) min{p, cr -1 }, c(f2) some positive constant depending only on the domain. 

c) The bilinear form b satisfies the inf-sup condition 

■ f b({v,K),(q,i>)) 

mf sup — -——- — - — - > 0. 

(v,if)ex (vK) ||(v,if)||x (ViKr) ||(9, V’)I|m ( 9iV , ) 
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Note that Hp(fl) (~l N(D) = {0}, for example, if Ti ^ 0. Otherwise, the subsequent analysis 
would require the introduction of a suitable quotient space of Hp(fl). 

Due to [7, Corollary 1] there exist the liftings of the boundary data 

u 0 e H 1 (f2) with V • u 0 = 0 and u 0 | Fl = gi, u 0 • n|p 3 = g 3 and u 0 - (u 0 • n)n|p 4 = g 4 , 

Jo € £ 2 (fl) with V • Jo = 0 and J 0 • n = j on E 4 . 

Define u := u u 0 , J := J — Jo- Also set <j>:= (j) — (f >o and p := p — p 0 , where po = 0, if r 2 UT 4 ^ 0, 
and po = Tprr / P, otherwise, and (f >o is the H 1 —lifting of the boundary data k on E 2 , if E 2 ^ 0, 

l“l Jn 

If. 

and (j )o = yp-y / (f >, otherwise. Define the bilinear form a : X( v x X ( v k) ► R and the linear form 

l“l Jn 

l '■ X(v,k-) ^ R by 

o((vi, Ki), (v 2 ,K 2 )) := a 0 ((vi,Xi), (v 2 ,K 2 )) + a 1 ((v 1 ,K 1 ), (uo, J 0 ), (v 2 ,K 2 )) 

+ai((u 0 , J 0 ), (vi, Jfi), (v 2 , K 2 )) 

and 

£(v, K) \= x(v, K) - a 0 ((u 0 , J 0 ), (v, K)) - ai((u 0 , J 0 ), (u 0 , J 0 ), (v, K)) - 6((v, K), (p 0 , </> 0 )). 
Note that if (uo, Jo) have small norm in X( u j), then the form a is coercive on X( V) jf), i.e. 

a((v,K),(-v,K))> .(a 0 -2||o 1 || ||(u 0 , Jo)||x (u>J) ) ||(v, (H) 

for all (v, K) £ X (v K y Substituting u = u + uo, J = J + Jo, (/> = (/> + 4>o and P — P + Po into 
Problem 1 we get its equivalent formulation. 

Problem 2: Find (u, J) € X( v x) an d (p, <j>) £ satisfying 

a{(u,J),(v,K)) +«i((u, J),(u, J),(v,K)) +b((v,K),(p,4>)) =£(v,K) 
for all (v,K) eX (V)K) and &((u, J), {q, if)) = 0 for all (q,ip) £ M^y 
Next, define 

V := {(v,K) £ X (V<K) : &((v, K), (q, if)) = 0 for all (qrf) £ M M) ] 


and consider 

Problem 3: Find (u, J) € V such that 

°((u, J), (v, K)) + ai((u, J), (u, J), (v, K)) =£(v,K) 


for all (v, K) £ V. 

The standard nonlinear version of the classical LBB (Ladyzhenskaya-Babushka-Brezzi) theory (see 
Chapter 4.1 in [26]) tells us that Problem 2 and Problem 3 are equivalent, if the form b satisfies the 
inf-sup condition, and that Problem 3 is uniquely solvable, if the form a is coercive and bounded. 
In other words, due to Lemma 3.1 and (11), the LBB-theory allows us to further transform Problem 
2 and solve Problem 3 for the unknown velocity u and electric current density J on a subspace V of 
X( v K y Thus, if (u, J) £ V is a solution of Problem 3, then there exist a unique pair (p, </>) € 
such that (u, J,p, iff) solves Problem 2. And, given any ((u, J), (p, <j>)) £ X( v ^ x M( q ^, solution 
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of Problem 2, then (u, J) is in V and solves Problem 3. By [7, Theorem in Section 3] Problem 
3 is uniquely solvable if 

||(uo, Jo)||x (UiJ) <^j (12) 

and 


W\W> < 


(« 0 - 2||« 1 ||||(uo, Jo)||x (UiJ) ) 2 
4||ai|| 


(13) 


The definition of ao = c(fi) minj^, a *} implies that any given set of data will satisfy (12)-(13) if 
the viscosity 77 and the electric resistivity a^ 1 are large enough. 


4 From weak formulation to frame discretization 

We start by reformulating Problem 3 to fit the following abstract setting: There exists a separable 
Hilbert space 7 ~L such that V C H C V' with bounded and dense inclusions. The triple (Y,H, V') 
is then called a Gelfand triple. The duality between V' and V is identified on H using the inner- 
product (-, -)u of H. There exists an operator A : V —> V' such that {Av,w)v x v '■= a(v,w) defines 
an elliptic bilinear form, i.e., there exist positive constants a, (3 such that a||u||y- < a(v,v) < 0\\v\\ v 
for all v G V. The cllipticity of a implies that ||Hu||v' ~ f|u||v and that A is a boundedly invertible 
operator with ||H _1 1| < a -1 . We also assume that there exists a trilinear form ai inducing a bounded 
bilinear operator A\ : V x V —> V', defined by (Ni(i>, w), z)v xv := ai(v,w,z) for v,w,z £ V. It 
has been shown in [7] that Problem 3 translates into the following abstract problem: 

Find u := (u, J) £ V such that 


Au + A\ ( u, u ) = (14) 

where l £ V' is a functional on V. 

The solvability of (14) is ensured by [7, Lemma 2] and summarized by the following theorem. 
Theorem 4.1. Let H : V —> V be given by H(v) := A -1 (£ — Ai(v,v)). The following hold 

(a) If 0 < r < 2 \\a 1 \\ ’ then H\b t is a contraction with Lipschitz constant L := 2 ^ 1 ^ r , where 
B r C V is the closed ball of radius r centered at the origin; 

(b) I/O < r < -pYjj and ||f|| v' < r(a — ||^4i||r*), then H\ Br : B r — » B r ; 

(c) If ||^||v' < 4 p 7 [f’ then (14) has a unique solution u with ||it||v < 2 lfT 7 Tf an< ^ so tufion is 

given by the fixed point iteration 

u n+ 1 = Hu n , u 0 = 0, n G N 0 , (15) 

u = lim u n . 


Note that (15) is equivalent to 


Au n+ i=£-A 1 (u n ,u n ), n G N 0 , u 0 = 0. (16) 

Under our assumptions on A , the equations in (16) are elliptic operator equations. We show in 
this section how to derive a discrete problem equivalent to the abstract nonlinear problem in (14) 
and prove a result similar to Theorem 4.1 showing that the solution of the discrete problem exists 
and is unique under certain assumptions on the parameters of the original MHD problem. The 



discrete problem is obtained using suitable stable, redundant, and nonorthogonal expansions, so- 
called Gelfand frames for the Gelfand triple (V, TL, V'). We also show that the corresponding discrete 
fixed point iteration can be numerically realized efficiently. In particular, we present the realization 
of the key numerical routine SOLVE used in our discrete fixed point iteration to approximate 
adaptively the solution of the elliptic problems in (16). 

4.1 Gelfand Frames 

In the following, the sequence space £2 (A/ - ) on the countable index set Af C R d is induced by the 
norm 

l c llr 2 (A0 • ( ^ , |Cn| J ; c = {Cri}neA/" € £2{Af). 

VnG Af ) 

The space £q (Af) C £2(Af) is the subspace of sequences with compact support. Denote (•> •)n and 
|| • || the inner product and the norm on the separable Hilbert space TL, respectively. A sequence 
£F := {/njngy in 'H is a frame for TL if 

ll/ll^£|</,/nM 2 , for all fen. (17) 

neAf 

Due to (17) the corresponding operators of analysis and synthesis given by 

F:H^£ 2 ( AT), f~({fJn)n) neAr , (18) 

F* : £ 2 {Af) —> TL, c^£c„/ n , (19) 

are bounded. The composition S := F*F is a boundedly invertible (positive and self-adjoint) 
operator called the frame operator and T := {S^ 1 f n }n&Af is again a frame for TL, called the canonical 
dual frame, with corresponding analysis and synthesis operators 

F := F(F*F)~ 1 , F* := {F*F)~ 1 F*. (20) 

In particular, one has the following orthogonal decomposition of £2{Af) 

£ 2 (AT) = ran(F) ® ker(F*), (21) 

and 

Q := F(F*F)~ 1 F* : £ 2 (Af) -f ran(F), (22) 

is the orthogonal projection onto ran(F). 

The frame T is a Riesz basis for TL if and only if ker(F*) = {0}. In general, we assume that 
{0} is a proper subspace of ker(F*). In other words, due to the redundancy of the frame there 
may exist sequences c = {c„} n e jv d = {ci„}neA r hi £2(Af) such that £ Cnfn — £ dn fn • la 

n£jV neAf 

particular, the redundancy may lead to the situation when a small perturbation d of the coefficient 
sequence c has no effect on the synthesis operator. This possible reduction effect on errors, noise, 
and numerical round-offs is the motivation for using frames for the applications, where tolerance to 
errors is required. The intrinsic stability of frames is also expected to play a role in the conditioning 
of the discretizations of operator equations and leads to additional robustness that the discretization 
inherits from the frame. In fact, it has been recently confirmed by numerical experiments in [48] that 
increasing (uniform) redundancy of the frame one improves the conditioning of the corresponding 
discretization matrices. 
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It is somewhat perplexing that different sets of coefficients yield equivalent representations of 
the same element of TL. It is not clear then what are the “good and computable coefficients”: The 
importance of the canonical dual frame is its use in reconstruction of any f £ Ti , i.e. 

/ = SS-\f = J2 if, S^fnUfn = S~ 1 Sf = J2 </, fnUS-'U (23) 

n€Af n£Af 

Since a frame is typically overcomplete in the sense that the coefficient functionals c = {c n (f)} n ^j\f £ 
£2 (A/”) in the representation 

/ = E c « (/)•/» ( 24 ) 

neAf 

are in general not unique (ker(F*) ^ {0}), there exist many possible non-canonical duals {f n }neAt 
in 7 i, for which 

(25) 

nGAf 


A more general definition of frames is required for Banach spaces. For details on Banach frames 
we refer, for example, to [13, 24, 25, 28, 30]. Throughout this paper we make use of Gelfand frames 
that are particular instances of Banach frames. So we do not introduce the latter here in full 
generality. Assuming that B is a Banach space continuously and densely embedded in 'H we get 

B c H ~ n' C B’. (26) 

If the right inclusion is dense, then (B, H,B') is called a Gelfand triple. The symbol ~ stands for 
the canonical Riesz identification of 'H with its dual TC. 

Deflnition 1. A frame T (here T is some dual frame, e.g., the canonical dual frame) for H is called 
a Gelfand frame for the Gelfand triple (B, H, £>'), if T C B, T C B’ and there exists a Gelfand triple 
(Bdi £ 2 (A/"), B ' d ) of sequence spaces such that 

F*:B d B, F*c= ^ c n f n and F : B -> B d , Ff = ((fJn) B xB') neM (27) 


are bounded operators. 


REMARKS: 

1. If T (again T is some dual frame, e.g., the canonical dual frame) is a Gelfand frame for the 
Gelfand triple ( B,H,B ') with respect to the Gelfand triple of sequences (B d , ^(A/”), B ' d ), then by 
duality also the operators 

F* : B' d - B\ F*c = ]T c n f n and F : B' —> B' d , Ff = ((f,f n } B > xB ) n ^ (28) 

nGAf 

are bounded, see, e.g., [31] for details. 

2. If B = R then Definition 1 becomes the definition of frames for Hilbert spaces. 

3. Even if the solution space V C TL = L 2 (fl) C V' is a Hilbert space, by using the notation 
U B" here we want to emphasize that the frame T we consider is not a Hilbert space frame for V. 
It is a frame for 7f, which also characterizes V (as a subspace of 1~L) with frame coefficients (/, fi)u 
belonging to some suitable sequence space C £2 (A f) and being computed using H -inner product. 
Of course, we can consider genuine Hilbert frame expansions for V (as it is done, for example, in 
[40]). This may create some confusion because one might be tempted to use the V—inner product 
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when determining the duality pair. This, in most of the cases, is not compatible with numerical 
implementations. 

4. Definition 1 generalizes the following case to pure frames: Consider a wavelet system 4' := 
{4>j, k}j>-i,kejj on n (Jj is a suitable set of indexes depending on the scale j, see [46] for details), 
B = H s (nj,H = L 2 {ty and 


Bd = £ 2 , 2 s := {d {dj,k}j>-i,kejj 



< 00 }. 


It is well known that if 4/ is a Riesz basis for and its elements, together with those of its 

biorthogonal dual basis 4/ := > are compactly supported, smooth enough, and with 

a sufficient number of vanishing moments, then is fully characterized by 41 in the sense that 

/ € H s ( f2) if and only if 

f = El (f L 2 (Cl)' t Pj,k (29) 

j : > 1 fcey,’ 


and 


1/2 


I H*(Q) 


E E 22 sj K/,E 


k/L 2 (n) I 


US —1 k£jj 


(30) 


See [13] for the same characterization by using pure wavelet frames, constructed by Overlapping 
Domain Decomposition. Note that there exists a natural unitary isomorphism from £ 2 , 2 s - into £2 
given by 


D 




■ £2 


£ 2 , d := {dj,k}j>-\,kejj -Dff s (f2)d := {2 Js dj^}j>-i,keJj- 


(31) 


Keeping in mind the example in the above REMARK 4., we proceed to the numerical treatment 
of abstract elliptic operator equations by means of Gelfand frame discretizations. 

4.2 Adaptive Numerical Frame Schemes for Elliptic Operator Equations 

To implement the fixed point iteration described in Theorem 4.1, we first study the solvability (for 
fixed u^"i) of the linear operator equations in (16). Generally, such equations are of the form 

Au = /, (32) 

where A, as before, is a boundedly invertible operator from Hilbert space V into its dual V', 

||Aw||v' ~ IMIv, u e V. (33) 

We also have that 

a(v,w) := (Av,w) V 'xv , (34) 

defines a bilinear form on V, where (•, -)vxV' defines the dual pairing of V and V'. The form a is 
elliptic, i.e., there exist positive constants a, (3 such that 

a IMIv — a i v > v ) < /?IM!v (35) 

for all v € V, and a is non-symmetric. The assumption on the non-symmetry of a is motivated by 
the MHD example presented in Sections 2-3. 

Here and throughout the rest of the paper we assume that T = {/ m } n ey is a Gelfand frame 
for the Gelfand triple (V, W,V') with (Vd,£ 2 {Af),'V' d ) being the corresponding Gelfand triple of 
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sequence spaces. Moreover, keeping in mind REMARK 4. in Subsection 4.1, we also assume that 
there exists a unitary isomorphism TV : Vd —> t? 2 (A/"), so that its £ 2 (A/")-adjoint D^ : Z 2 (A/”) —> 
is also an isomorphism. 

We show, next, how the Gelfand frame setting can be used for the adaptive numerical treatment 
of elliptic operator equations (32). Following, e.g. [8, 13, 40], one uses frame expansions to convert 
the problem (32) into an operator equation on £ 2 (A/"). The problem that arises is that the redundancy 
of the frame leads to a singular discretization matrix. Nevertheless, in Theorem 4.3 below we show 
that this can be handled in practice and that the solution of (32) can be computed by a version 
of Richardson iteration applied to the associated normal equations. The resulting scheme is not 
directly implementable since one has to deal with infinite matrices and vectors. Therefore, similarly 
to [8, 9, 40], we also show how the scheme can be transformed into an implementable scheme using 
“finite” versions of the building blocks procedures we introduce in Subsection 4.2.2. The result is a 
convergent adaptive frame algorithm. 

4.2.1 A series representation 

We start by generalizing Lemma 4.1 and Theorem 4.2 given in [13] to the case of non-symmetric a. 
We give the detailed proofs to emphasize the difference between the symmetric and non-symmetric 
cases. 

Lemma 4.2. Under the assumptions (34), (35) on A, the operator 

A := (D^)~ 1 FAF*D~ 1 (36) 

is a bounded operator from £2 (Af) to £ 2 {Af). Moreover A is boundedly invertible on its range 
ran(A) = ran((D(»- 1 F). 

Proof. Since A is a composition of bounded operators : ^ 2 (A/") —> Vd, F* : Vd —> V, A : V —> 
V', F : V' —> and (Z?^) -1 : —> ^ 2 (A/ - ), A is a bounded operator from £ 2 (Af) to £2 (A/”). 

Moreover, from the decomposition (36) we get 

ker(A) = ker (F*D^ 1 ), ran(A) = ran((D^) _1 F). (37) 

Define L := F F* Dfj 1 . Note that ker(L) = ker '(F*D^ 1 ) and ran(L) = ran((D^) _1 i 7 ’). The 

fact that ^ 2 (A/") = ker(L*) © ran(L) implies, due to the self-adjointness of L, that 

e 2 {AO = ker (F*D^) © ran((D^)- 1 F). (38) 

Therefore, 

Aj ran(A) : ran(A) -> ran(A) (39) 

is boundedly invertible. B 


Denote by P : £ 2 (A/") —> ran(A) the orthogonal projection of £ 2 (A/") onto ran(A). 
Theorem 4.3. Let A satisfy (34) and (35). Denote 

f :={D* w )~'Ff 

and A as in (36). Then the solution u of (32) can be computed by 


u = F*D w l Pu 


where u solves 


Pu = (ex* J2 (id —a*A*A)" ran(A) A*f, 


n—0 


with 0 < a* < 2/A max , where A ma x = || A*A ||2 with || • H 2 being the usual spectral norm. 


(40) 

(41) 

(42) 
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Proof. We have u = J2neM^ u ' fn)nfn in H. Since T is a Gelfand frame, F*F : V —> V is bounded 
and implies u = F*Fu = YlneAf( u i /JvxV'/n in V. Moreover, (32) is equivalent to the following 
system of equations 

{ u i /n)vxV' (Af n , /m)y'xV = (/i/ro)v'xVi m € Af. (43) 

neAT 


Denote u := D^Fu and f, A as in (40) and (36). Then (43) can be rewritten as 

Au = f. (44) 

Multiplying both side of (44) by A* we get the normal equation 

(A*A) u = A*f. (45) 


Note that A* A is self-adjoint and positive-definite by the hypothesis. Note also that ker(A*) is 
orthogonal to ran(A) and (38) implies that ker(A*) = ker(A). This and the invertibility of A on 
its range implies that A*A is boundedly invertible on ran(A). Therefore, the solution of (44) is 
equivalent to the solution of (45). 

Note that, for 0 < a* < 2/A max , the operator 


B:=a*£(id-a*A*A)f ran(A) . 

n=0 


(46) 


is well-defined and bounded on ran(A), since p(a*) := || (id — a*A*A) ran ( A ^ 11 2 = max{a*A max — 
1,1 —a* A m in} < 1, where A m i n := || (A*A| ran ( A )) _1 1 |2 ■ The function p is minimal at a* := 2/(A ma x + 
Amin)- Moreover, 

B O (A*A| r an(A)) = (A*A) O B| ran (A) = id| ran(A) • (47) 


Since A(id—P) = 0, 


Au = APu = f. 


(48) 


Therefore Pu G ran(A) is the unique solution of (44) in ran(A) and by (47) 


By construction 


{f, fm)v'xv 


so that u = F*D V 1 Pu solves (32). 


Pu = BA*f. 


( F*Ff , / m )v'xv 
{F*D^f, / m ) V 'xv 
(F*DvAPu, /rn) V 'xV 

(AF*D~ 1 Pu, fmW'xv, toga/ - , 


(49) 


4.2.2 Numerical realization 

Now we turn to the numerical treatment of (44). Due to Theorem 4.3, the computation of u solving 
(44) amounts to an application of the following damped Richardson iteration 

u(* +1 ) = u (i) - a*A*(Au (i) - f), i G N 0 , u (0) = 0. (50) 

Certainly this iteration cannot be practically realized for infinite vectors u^\ i G No- To avoid this 
problem, we make use of the following procedures (see [8, 9, 10, 40] for details on their analysis and 
numerical realization) : 


13 



(51) 


• RHS[e, f] -> f e : determines for f G ^(W) a vector f e G such that 

Ilf — ?e|U 2 (AG) < e; 

• APPLY [e, A,v] — > w e : determines for a bounded linear operator A on (^{M) and for v G 
to(N) a vector w e G to (AO such that 

IIAv — Well^^) <e; (52) 

• COARSE[e,v] — > v e : determines for v G ta{M) a vector v e G £o(Af) such that 

l|v- v e \\ t2{M) < e. (53) 

We discuss in more details further properties of the routines RHS, APPLY and COARSE in 
Section 6, where we study the complexity and the computational cost required to approximate the 
solution of the original problem up to some prescribed tolerance. 

Let p := p(a*) be as in the proof of Theorem 4.3. We can define the following inexact version of 
the damped Richardson iteration (50): 

Algorithm 1. 



The parameter 9 plays an important role in complexity estimates (given in Section 6) for 

COARSE. 

The proof of the convergence of Algorithm 1 is analogous to that of [40, Proposition 2.1] and of 
[13, Theorem 4.2] except for the fact that here we make use of the damped Richardson iteration in 
(50) on normal equations, due to the non-symmetry of a. Nevertheless, the following result holds. 

Theorem 4.4. Under assumptions of Theorem 4-3, let u G be a solution of (44). Then 

SOLVE[e, A, f] produces finitely supported vectors v^’K ), ^(i) ; such that 

||P(u - v «)|| £2(a0 < e,', j G No. (54) 
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Figure 1: Procedure SOLVE applied to the Poisson equation on the L-shaped domain to ap¬ 
proximate the solution u(r,9) := £(r)r 2 / 3 sin^#) + u(r,9), C G C 0O (fl) is a suitable truncation 
function corresponding to the reentrant corner of the L-shaped domain, and u £ H 2 (Q) fl 
Overlapping rectangular patches and aggregate wavelet frames are used for the discretization. See 
[13, 14, 40, 48] for more details on the numerical implementation. Approximations are shown for 
the case of piecewise linear wavelet frame elements. 


In particular, 


Moreover, it holds that 


\u - F*DF v, 


v 1 v e || v <||Fl|| J Dv 1 | 


||Pu - (id -P)^- 1 ) - v(^) || MA0 < -J-, j> 1. (56) 

Of course, the numerical implementation of the damped Richardson iteration on the normal 
equations (45) might exhibit a low convergence rate if the relaxation parameter a* is small. To 
improve the efficiency of the proposed scheme, the generalizations of Algorithm 1 towards, e.g., 
(conjugate) gradient iterations as suggested in [19, 48] are now a matter of investigation, see also 

[14]- 

5 Numerical realization of the fixed point iteration 

Now we have at hand the major building blocks needed to formulate an implementable fixed point 
iteration. In this section we show how the problem in (14) can be discretized and how Algorithm 1 
can be used to implement the fixed point iteration in (15). 

We denote 

i) A := (D^)^ 1 FAF*D^ 1 : £z(N) —> £ 2 (A/ - ), bounded and boundedly invertible on its range 
ran(A) = ran((DL) _1 F) ; 




























Figure 2: Reduction of the residual ^ 2 -norm and the L 2 -evror for SOLVE (dotted line) and for 
the adaptive steepest descent (solid line). We refer to [14] for more detailed numerical tests where 
optimality is shown and discussed. 


ii) 1 := G ran(A) C I 2 (N)\ 

hi) A,(*) := (D^)- 1 FA 1 (F*D~ 1 -,F*D V 1 -) : t 2 (Af) - ran(A) C £ 2 {M); 

iv) H(-) := (Ai^a)- 1 )?- Ai(-)) : £ 2 (Af) -+ ran(A). 

REMARK: Note that ran(A) = ran((£>^) _1 i 7 ') implies that the operators in ii)—iv) map I 2 (M) into 
ran(A). By definition of P we also have H(Pv) = PH(v) = H(v) for any v G £ 2 (J\f). 

Then, it is easily verified that the variational problem in (14) is equivalent to the following 
discrete problem. 

Problem 4: Find u G ran(A) C I 2 (A/") such that 

Au + Ai(u)=T, (57) 

or, equivalently, such that u is a fixed point of H in ran(A), i.e 

u = H(u), u G ran(A). (58) 


Define the closed subset of ran(A) by 

B r \= {u G ran(A) : ]|u||^ 2 (a/‘) < r} for some r G R + . 

The next theorem is a discrete analogue of Lemma 2 and Corollary 2 in [7]. 

Theorem 5.1. Under the assumptions and notations specified above, the following statements hold 
true: 

a) If 0 < r < (2|| (A| ran a) -1 || ||E|| 3 || ||) , then H| Br is a contraction with Lipschitz constant 

L:=r(2|KA jranA )- 1 ||||E|| 3 ||A 1 ||)<l; 

b) if 0 < r < (|| (A| ran a) 111| A 1 II 3 II Ai ||) 1 and ||1|| < r (||(A| ranA ) - 1 |r 1 -|^i||||^l| 3r ) then 

H (B r ) C B r ; 

c) if ||1*|| < (4||(A| ran A)^ 1 || 2 ||T'|| 3 ||Ai||) 1 then (58) has a unique solution u G B r *, for some 
suitable r* such that 0 < r* < (21| (A| ran a) — 1 IIII11II II) 
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Proof. Recall that Z?v is assumed unitary and ||-Dv|| = H-Dv 1 II = ll^vll = 1- Moreover, since F 
maps V' into both being Hilbert spaces, thus ||F|| = ||F*||. Then, for u, v e B r , 

l|Hu — Hv|| < ||(A |ranA )- 1 ||||A 1 u - 

= II (A, ran a)' 1 IIII (D^-'F (A^D^u, F*D^ u) - A\(F*D^v, F*D^v)) \\ UM) 

< II (A, ran a)' 1 II l|F|| || A^D-'u, F*D~ 1 u) - A 1 (F*D~ 1 v, F*D^v) || V ' 

= l|(A| ranA )- 1 ||||i ;i ||p 1 (F*^- 1 (u-v),F* J D- 1 u)-H 1 (F* J D- 1 v,F^- 1 (v-u))|| v , 

< ||(A| ran A)- 1 ||||^||||^ 1 ||||^|| 2 ||i3 -v|U 2(AA) (||u|U 2( aA) + l|v|U 2(J v)) 

< 2r||(A| ran A) _1 ||||F||^|Hi||||u- v|| MA0 . 

Thus, as L < 1 we have that H is a contraction. To show b), just observe that by definition of H 
and the estimate above 

||Hu||, 2(a0 < ||A^ anA || (ill'll + ||H 1 ||||F||V) <r. 

c) We have to show that if ||1|| < (4||(A| ran A) _1 || 2 ||F|| 3 ||A 1 ||) 1 then there exists r* with 0 < r* < 

(2||(A| ran A) _1 ||||F|| 3 ||A 1 ||) 1 such that ||f|| < r* (||( A, ranA ) _1 II~ 1 - II Ai || ||F|| 3 r*) . Then, using 
parts a) and b) of this Lemma we get that H\b t * is a contractive mapping of B r * into itself and, 
thus, has a unique fixed point. To see that such an r* exists, consider 

h(r) = r (||( A, rarL A ) 1 II 1 - II Ml || ||F|| 3 r) , 

a quadratic mapping, and note that h assumes values from 0 to (4||(A| ranA ) _1 || 2 ||F|| 3 ||A 1 ||) 1 as r 
varies from 0 to (2||(A| ran A) _1 |]|F|| 3 ||Ai||) \ ■ 

Corollary 5.2. If ||f|| < (4||(A| ran A) -1 || 2 ||F|| 3 ||Ai||) -1 , then the solution u in B r * of (58) is given 
by the following discrete fixed point iteration 

fi„ + i = Hu„, n > 0, flo = 0 S ran(A), (59) 

u = lim u„. (60) 

n —»oo 

The discrete fixed point iteration cannot be implemented for infinite vectors. For this reason we 
propose the following new approximating adaptive scheme FIXPT, where the procedure SOLVE 
introduced in Subsection 4.2.2 replaces the exact computation of u n+ i in (59). 

Algorithm 2. 
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REMARK: The Algorithm 2 is a perturbation of the exact fixed point iteration creating sequences 
{vj}ieN 0 of finite vectors. In general such vectors will not belong to ran(A) and there is not much 
hope that limn-,,^ v n = u. However, one might try to see whether limn-,^ Pv„ = u. A priori it 
might even happen that, because of an accumulation of the perturbation errors, there exists some 
n G N large enough such that Pv„ ^ B r • ! In order to show the convergence of Algorithm 2 to the 
solution of problem (58) we prove the following Lemma. 

Lemma 5.3. If ||1|| < (4||(A| ran A) _1 || 2 ||-P|| 3 ||Ai||) _1 and if the quantities 

n+1 n—3 n—h n— 1 

£n ■— + (1 + L) £kL n h k + ^£l + L(so + || (A| ran (A)) 1 II l|lj|)') L k 

fc=2 h =0 fc=3 fc=0 

+ e 0 + ||A^ (A) ||||Ti| < r*, for all n G N 0 , 

with L and r* as in Theorem 5.1 a) and c), respectively, then the sequence {Pv+eN resulting from 
the application of Algorithm 2 all lies in B r .. 

Proof. We show by induction on n that 


||P(v n +i — v„)lk(A 0 < £ n+i + (1 + L) k + L n 1 (ei + L||Pvi||£ 2 (a^)). (61) 

k —3 


and 


n+1 


n—3 n—h 


2—h—k 


|P(v„+i)lk(A 0 Y £k + (1 + L) Y. Y, £ kL 

k—2 h—0 k =3 

+ (ei + L(e 0 + ||(Aj ran ( A )) 1 ||||1||)) Y,L k + £ o + l|A| ra 


(62) 


n—1 


ran(A) I 


k =0 


Set n = 1. We get 

l|Pvi||^(^) = ll p (sOLVE[ £o ,A,i]-H(0))+H(0)||, 2(A , ) < £ 0 + ||A 


- 1 ] 
| ran(A) 


< £ 0 + II (Aj ran (A)) 1 1|||1|| < -r*. 


Moreover, since H(v) = H(Pv) = PH(v), we get 


P (SOLVER, A.l-A^Vi)]) - H(v!) = P (SOLVER, A.l-A^vl)]) - P H( 


vi). 


By (54) it holds that ||P SOLVE[£ n+ i, A, 1—Ai(v n )] —P H(v„)|| < £ n+ i and, using the assumption 
on ||1|| ensuring that H is a contraction with the Lipschitz constant L , we get 


||P(v 2 -Vl)|.k(A 0 


= P (sOLVE[£ 2 , A,f-Ai(vi)]) -P (sOLVE[£i, A,f-A : 
= S p (s°LVE[£ 2 ,A,r-A 1 (v 1 )] -H(vi)) +H(v 1 )-H(v 0 ) 
+ H(v 0 )-P (SOLVER, A, t-Ai+o)]) Ik (AO 

< £ 2 + £i + L||Pvi|| f2 (A/-). 


4>(A 0 
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We assume now that all Pvi,...,Pv n £ B r * and that formulas (61) and (62) are valid replacing n 
with n — 1. Then 


l|P(Vn+l - V n )|k(AT) < £n+l+e n +L\\P{^n-^ n -l)\\e 2 (M) 

( n—1 

£ n + (1 + L) £kL n 1 k + L n 2 (ei + LllPV! ||£ 2 (A/')) 
fc=3 
n 

= £n+l + (1 + i) £kL n k + L n 1 (£l + T||Pvi||^ 2 (^)). 

k —3 

Using the above estimate for ||P(v n+1 — v+I^aT)) the triangular inequality, i.e. 

ll P Vn+l||<! 2 (A0 ^ || P (V „+1 - V n )|| fe(A/ -) + ll P + \\t 2 (AT) 

and (62) for Pv„ we have immediately (62) for Pv n+1 . Thus, by hypothesis on £ n we obtain 
|jPv„ +1 ||^ 2 ( 7 y) < r*. This implies by induction that Pv„ e B r . for all n € N. I 


REMARK: Note that the assumption on £ n in the above Theorem is not restrictive as soon as the 
data are small. In fact, it holds that 


n+1 

E 

k—2 


n+1 

£ k = E 

k =0 


£q ~ 1 - £ 0 = 


I _ £ n + 2 
1 b 0 

1 - £o 


- 1 - £0 < 


1 — £q 


0 for £q —> 0; 


• the map 

n—3n—h 

h —0 k—3 

where C(£q,L) —> 0 for £q —► 0, uniformly with respect to n; 


• (ei + L(£o + II A! r ^ n(A) || ||f||) 5E + (£o + II A! r ^ n(A) II ||T||) is small as one wants whenever £ 0 

fc=o 

and ||1|| are small enough. Note that ||1|| < ||F||||^||, where ||£||, at least for the MHD problem 
discussed above, depends on the size of the norms of the forcing terms and boundary data 
in (l)-(9) by means of formula (13). Unfortunately for nonlinear problems in fluid dynamics 
existence and uniqueness of the solutions are ensured only under smallness assumptions as 
in (13), so that we cannot expect that numerical implementations can be less restrictive. 
Nevertheless, such restrictions represent as usual the worst case analysis. 


We are now ready to show our main convergence result. 

Theorem 5.4. If ||1|| < (4||(A| ran A) -1 || 2 ||-F||' 5 ||+||) -1 and £„ < r* for all n € N, then 

H(u) = u = lim Pv„ = P (FIXPT[0, A, Ai, f]') , (63) 

n—>oo V / 


where {v.+ew are obtained by Algorithm 2. After n iterations of Algorithm 2 one gets 

*S<*o-l(+) +»-+)”) 

|Pv„+i - u|| WJ V) < £„-;- + £"INI < C»- - + LV. 


£o — L 


£o ~ L 


Therefore, for e > 0 one has 


||u-P(FIXPT[£,A,Ai,f]) || <£, 


(64) 


(65) 
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and the number N of iterations to achieve the accuracy e > 0 is estimated by 


AT ~ — log(e). (66) 

Proof. First we want to show that {Pv+ e N is a Cauchy sequence in B r *. To do that, we use the 
estimation (61) and get 


n+m m n-\-m—h 

||P(vn +m -v „)||, 2(A0 < ]T Sk + il + L)^ E L n+m ~ h ~ k 

k—n -\-1 h— 1 k—3 

n+m —2 

+ + L(eo + ||(A| ran (A)) Ullll)) L k . 

k—n—1 

A straight forward computation yields that the expression on the right in the estimate above goes 
to zero as n goes to infinity. Therefore, for any e > 0 there exists n £ N large enough such that 
for all to € N, to > n it holds that ||P(v ra+m — v„)||^ 2 (A0 < £ ■ Due to Lemma 5.3 the sequence 
{Pv,;}, 6 n £ B r *. Therefore, {Pv+gN is a Cauchy sequence in B r » and then (because B r » is a 
closed subset) there exists a unique v £ B r * such that v = lim Pv„, Moreover, we have 


Pvn+i = P (SOLVE[e„ + i, A, 1 - Ai(v„ 

= (p(sOLVE[ £ „ +1 ,A,r—Ar(v n )]) - H(Pv n )) + (H(Pv n ) - H(u)) + u. 


This implies that 


|Pv„+i -u|| = ||p(sOLVE[e n+ 1 ,A ) l-A 1 (^ n )]J - H(Pv n )||, 2w + ||H(Pv n ) - H(u)||, 2(AA) 

< £ n +l + L||Pv„ - u||^ 2 (Af) < £n+ 1 + L(e n + L||Pv„_i - u||<? 2 (a/-)) 

n 

£ o^£o fe L fc + L n 11u|U 2 (AT) 


< 


= £o 


= £o 


fe=o 

( e o +1 — L n+1 ) 


£o — L 

( / \ n \ 


+ L"||u| 


<b(A0 


l "I|u|+aO -»• 0, n —> oo. 

(' 


so — L 

Therefore, one has H(u) = u = lim^oo Pv ra = P ^FIXPT[0, A, Ai,lj) . Moreover, since 11u|U 2 (a^) 
r*, then the second inequality of (64) is valid and one has 


< 


£ o 


/ \ n 

£ o( £ o-L(±) ) 


if and only if 


£n 


£0 — L 
£q — L 


+«-+)) 


+ L n r* < £ 


C e-L n r *), 


(67) 


( 68 ) 


which is the criterion to exit the loop in Algorithm 2. This implies immediately (65). From (67) 
one shows that the necessary number of iterations of achieve accuracy £ > 0 is given by (66). H 


Similarly to the proof of formula (55) (see [13, Theorem 4.2]) one finally shows the following. 
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(69) 


Corollary 5.5. If ||1|| < (4|| (A| ranA ) 1 || 2 ||-F|| 3 ||Ai||) 1 and £ n < r* for all n G N ; then 

« = E (£>v lFIXp T[0»A ) A 1 ,i]) n / n , 

n£Af n 

is a solution in V of the fixed, point problem H ( u ) = u and one has that for e > 0 

E (r>v lFIXp T[e.A,A 1 ,i]) n / n ||v<e. 


u - 


r 


Proof. Since ’kei(F*D^ 1 ) = ker(A) = ker(P) and F*F = id, we finally verify 

\\u-F*Df / 1 u £ \\ v = ||P' t (Fw-^v lfl e)||v 

= ||P^v 1 (Pfl-fle)|| v 

= ||F^v lp («-u £ )|| v 

< ||F*||||Dv 1 ||||P(u-u £ )||, 2W . 

The claim follows from the above Theorem. 


(70) 


REMARK: The fact that at each iteration we compute the approximation up to a perturba¬ 
tion/tolerance £j also means that the scheme is stable. In other words, not only can £, be interpreted 
as the numerical approximation accuracy we achieve at each step, but also as the error tolerance we 
can afford, without spoiling convergence. Moreover, the scheme is fidly adaptive in the sense that 
the iterations are enforced (by the use of suitable implementation of COARSE, see below) to work 
only with minimal number of relevant quantities (frame coefficients), in order to keep the prescribed 
accuracy-complexity balance. 


6 Quasi—optimal complexity of the algorithm 

In this section we present the complexity analysis of Algorithm 2. From Theorem 5.4, in particular 
from formula (66), we already know that to achieve a prescribed accuracy e > 0 one needs to execute 
N ~ — log(s) iterations. Therefore, having an estimation of the cost of each iteration, the asymptotic 
analysis of the complexity of the suggested adaptive scheme can be done. This is one of the very 
interesting theoretical advantages of the adaptive (wavelet) frame approach, together with the fact 
that one can prove both convergence and stability of the adaptive scheme as shown in the previous 
section. 

Of course, the main ingredient of iterations in Algorithm 2 is the procedure SOLVE. As an¬ 
nounced in Subsection 4.2.2, we discuss here an implementation of such a procedure and study its 
complexity. To this end, one should illustrate how the building block procedures RHS, APPLY, 
and COARSE can be implemented and estimate their computational cost. Therefore, in this sec¬ 
tion we focus on the main properties and requirements of these building blocks so that Algorithms 1 
and 2 have certain complexity. We refer the interested reader to [9, 10, 13, 40] for the descriptions of 
procedures RHS, APPLY, and COARSE that fulfill the requirements stated below and needed 
for the complexity estimates. The complexity estimates for even more general algorithms than Al¬ 
gorithm 2 for linear and nonlinear variational problems, but under the more restrictive assumption 
that the discretizing frame T is a Riesz basis, have been given in [9, 10]. In order to describe the 
complexity in the more general case of pure frame discretizations a bit more (technical) effort and 
preparation is needed. We start by defining a so-called sparseness class of vectors, s € R + . A s 
will turn out to be such that, if u £ A s , then the size of the support of u e = FIXPT[e, A, Ai, l] 
and the computational cost for obtaining u £ can be estimated a priori. 
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An algorithm for computing a finite approximation u £ of u up to e, for u given implicitly as a 
solution of some equation, is called optimal if # supp(u e ) (the number of elements of the support of 
u e ) is not asymptotically larger (for £ —> 0) than the same quantity obtained by direct computation 
of any other approximation of u using the same tolerance £, for u being given explicitly. In addition 
to this, the optimality is fully realized if the complexity to compute u e does not exceed # supp(u e ) 
asymptotically (for e —> 0). In other words, one does want the number of algebraic operations to 
be comparable to the size of what is being computed. For discretizations by means of Riesz bases 
optimal algorithms can be realized, see [8, 9, 10], for example. Analogous algorithms for frames 
may exhibit “arbitrarily small reductions” with respect to the expected optimality (see the following 
Theorem 6.1 for the precise statement). Since the techniques for estimating complexity we use are 
similar to those in [40, Theorem 3.12], we encounter similar difficulties to achieve the full optimality 
for FIXPT when dealing with pure frames. We call this situation quasi-optimal. 

An optimal sparseness class is modeled as follows: For given s > 0 we define the space 

K,eak ■= {c G ^ 2 (AO : \\c\\ Ataeak ■= sup n 1/2+s | 7 „(c)| < oo}, (71) 

raSN 

where 7n(c) is the n— th largest coefficient in modulus of c. It turns out that |j • ||^t« is a quasi- 
norm and we refer to [8, 21] for further details on the quasi-Banach spaces A s weak . Such spaces can 
be usually found in literature under the notation (weak-£ r ) where r = (1/2 + s) -1 G (0, 2), and 
they are nothing but particular instances of Lorentz sequence spaces. Let us only mention that 

l|c|Uj, eofc ~ supiW||c-cW|k ( AO, (72) 

iVeN 

where cn is the best N— term approximation of c, i.e., the subsequence of c consisting of the N 
largest coefficients in modulus of c. In particular, (72) implies that for all £ > 0 there exists N e > 0 
large enough such that for all c the best N e —term approximation cjv e has the following properties 

(i) ||c - Civj£ 2 (A 0 ^ £ ; 

(ii) #supp(cArJ < £ _1/s ||c||J/£ eafc ; 

(hi) l|civJUj, eofc < ||c|U^ afc . 

Furthermore, from (72) one gets the following useful technical estimate 

l|c|U Lofc < (#supp(c)) s -icm 1)eafe (73) 

whenever 0 < s < s and c has finite support. This optimal class of vectors perfectly fits with 
complexity estimates for adaptive schemes for elliptic linear equations [8, 9, 13, 40]. For more general 
nonlinear problems a “weaker” version of the sparseness class has been introduced and denoted by 
A s tree in [10]. 

It is not known whether there exist frames for which the solutions of generic nonlinear equations, 
particularly for MHD equations, can have frame coefficients belonging to the sparseness classes 
described above. Therefore, we discuss here the requirements, fulfilled by A s weak and A| ree , that 
a generic sparseness class should have to ensure quasi-optimality of our scheme. And, in case 
the solution belongs to any sparseness class with such properties, the algorithm will behave as 
ensured theoretically. The numerical tests in [1, 46] for turbulent flows motivate our assumption 
that the (wavelet) frame coefficients of the corresponding solutions do belong to some of these generic 
sparseness classes (i.e., only few significant wavelet coefficients can be expected to be relevant in 
the representation of the solution). Our conceptual approach is also motivated by the need to 
simplify the presentation of our complexity result without going into rather technical details of the 
properties of the particular instances A s weak and A| ree of A s . We refer the reader to [10] for more 
specific details. 
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For s > 0 and for a nondecreasing function T : N —> N such that N < T(N) we call any space 
A s a T-sparseness class if A s weak C A s C A s weak and A s C A s for all s > s and if for all u€d s and 
for all e > 0 there exists a finite vector u £ with the properties 

a) ||u — u e || < e; 

b) T(#supp(u £ ))< £ - 1 / S ||u||^ s ; 

c) llwlU* <||fi|U-. 

In particular we assume that there exists a constant Ci(s) such that ||u + v||_ 4 <! < C'i(s)(||u||^« + 
Hvll^s). Of course, A s weak itself is a sparseness class with T = I. Moreover, there exist other T- 
sparseness classes different from A s weak , for example, the class Af ree dehned in [10, Formula (6.7)], 
that also turns out to be relevant in our context. 

Note that, for s > s > s > 0, by the inclusions A s weak c cd s c A s weak and (73) we have 

l|c|U«‘ < ll 6 IUj ieafc < (# s upp(5)) s - s ||c|U^ ofe < (T(#supp(c))) s - s ||c|U s , (74) 

for all finite vectors c. 

Now we are ready to formulate our main conceptual requirements. For a fixed s > 0 
(Al) Let 9 < 1/3. We assume that for any e > 0, v G A s and any finitely supported w such that 

||v-w|| < 0e, 

for w* = COARSE[(l — 9)e, w] it holds that 

T(#supp(w*))< £ - 1 / S ||v||^, 


and 

||w*|U. < C' 2 (s)||v|U», 

for some costant C s (s) > 0. Moreover, we assume that the number of algebraic operation 
needed to compute w* := COARSE[e, w] for any finite vector w can be estimated by 
< T(#supp(w)) + o(e _1 / s ||w||y/). See [40, Proposition 3.2] and [10, Proposition 6.3] for 
examples of procedures COARSE with such properties. As it will be clear in the proof of 
Theorem 6.1 the use of the procedure COARSE is fundamental to ensure that the supports 
of the iterates generated by the algorithm can be controlled. Due to the redundancy of the 
system used for the discretization, one can expect that the introduction of greedy algorithms 
of matching pursuit type [20, 27, 33, 42, 43] can potentially allow even higher rate of compres¬ 
sion of the active coefficients with respect to the current implementations of COARSE. We 
postpone this investigation to forthcoming work. 

(A2) Thevectorw e := APPLYfe, A, v] is such that ||w £ ||^ s < ||v||^ s , T(#supp(w e )) < £ -1 / s ||v||y/, 
and it is computed with a number of algebraic operations estimable by < £ -1 / s ||v||y,f + 
T(#supp(v)). See [40, Proposition 3.8] and [10, Corollary 7.5] for examples of procedures 
APPLY with such properties. 

(A3) One of the most crucial procedures of the iterative approximate fixed point scheme is the re¬ 
alization of := RHS[ T |^,f- Ai(vj)] for each step i and j of the outer and inner loops, 
respectively. In particular, it requires computing efficiently a finite approximation to A-j (v,;), 
where A| is some nonlinear operator and v, is a given finite vector. In their pioneering work 
[10, 11], Cohen, Dahmen, and DeVore have found an effective way for solving this problem for 
multiscale and wavelet expansions. Note that the efficient evaluation of nonlinear functionals 


23 



(as the ones we consider for the MHD problem) on coefficients as in [10, 11] is valid as soon 
as the reference bases are multiscale, i.e., the size of the supports of the basis functions de¬ 
cays dyadically for increasing levels, and their corresponding coefficients can be structured in 
suitable trees. Moreover, one needs that such bases characterize Besov spaces by certain norm 
equivalences. In Section 7, we present a construction of suitable wavelet frames that enjoy 
the above mentioned properties and, therefore, allow for a straightforward application of the 
results in [10, 11]. 


We only mention Cohen, Dahmen, and DeVore results here and refer to [10, 11] for technical 
details (in particular see [10, Theorem 7.3 and Theorem 7.4]): there exists a procedure RHS 


such that |g^' ) - (1 - Ai(v, ; 


< 


9e 


T(#supp(gp ) )) 


< e~ 1/s \ 


,1/s 

l-Ao > 


|g. 


O’) i 


< 1 


» — 12 <xK' 

11 Vj 11 and the number of algebraic operations needed to compute gf 1 is bounded by < 
^ 1,S ¥i\\ 1 i: +T(#supp(vi)). 


In the following the subscript index i refers to the iterations in the outer loop of the fixed point 
iteration and the superscript j refers to the inner loop iterations in SOLVE. Moreover, e,; and €j 
refer to the outer and inner loop tolerances, respectively. All estimations below hold asymptotically 
for e —» 0 ( s as in Algorithm 2 ). 


Theorem 6.1. For 0 < s < s < s let A s be a T-sparseness class and u £ A s , the solution of (57) 
as in Corollary 5.2. Assume that 


(i) (A1)-(A3) hold for all s £ (0, s]; 

(ii) P is bounded on A* for all t £ (0, s]; 

(in) K > 0 and 0 < 9 < 1/3 in Algorithm 1 are chosen so that 


C' 1 (s)C' 2 (s)||id-P||(3p A 76') l/s " 1 < 1. 


(75) 


Here the norm |j id —P || is the norm of id —P as an operator on A s ; 
(iv) the constants L,£q and r* satisfy 


L < £o < r* 


and p 


-K 


L , 1 

7^L + S °- 2 - 


for some <5q > 0. 


(76) 


Then for any e > 0 and S > 0 such that s/s = 1 + 5, the finite vector u e := FIXPTfe, A, Ai, 1] 
satisfies 

a) Hu-Puell^) < e; 

b) #supp(u e ) < £-( 1 + l5 )/ s ||u||7 f<5 ^7 

c) the number of algebraic operations needed to compute u e is < e _ ^ 1+l5 ^' s ||u||^ l5 ^' s . 

Proof. For the proof of part a) see Theorem 5.4. Next, we show part b). Assume that {v.j}ieN 0 
is the sequence of vectors generated in Algorithm 2. We want to show that T(# supp(v.;)) < 
^ (1+5)/ 1u||^ +5)/s for i large enough. Since P is bounded on A* for all t £ (0,5], it is also 
bounded on A s . Therefore Pu £ M+ Then for ej > 0 there exists a finite vector (Pu) e . such 

that ||Pu - (Pu)^7 2(J v) < and #supp((Pu) ej ) < T(#supp((Pu) ej )) < eT 1/s ||Pu||^// < 
e” 1 ^||u|| 1 ( / <f- Therefore, by (74) we have 


s/s-i 


(Pu) e 


< 


U 




(Pu) e 


< 


U 


is/s-i. 


Pull 


< 


U 


is/s-i. 


u IU« = ll U 


,s/s 


(77) 
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for any s > s > s > 0. From (56) we get that 

llPvf - (id-P)v. 


0-i) ,^0.^)1 


< 


2 6>e,- 


u 2 (A') - 3 ( 78 ) 

with v\ x := H(vj_i). Due to (76) and (64), we obtain for any i large enough and some <5o > 0 that 

||Pu - Pvrik ( A 0 < l]H(u)-H(v i _ 1 )||, 2(A0 

< LUPu-Pvi-iD^^) 


< L 


£ o 1 ( £ ’o ~ L{L/e 0 ) 1 2 ) | ri _ 2 ^* 


£o — L 
L ( 4(e 0 - L(L/£ 0 y- 2 ) 


£o 


Eo — L 


+ L l ~ r* 


+ L l ~ 2 £ 0 r* 


= £ o- 
£ o 


< 4 


L ( (£o - L(L/£ 0 y- 2 ) ( L 


£q — L 


£q — L 

f S 0 


V £ 0 


1 — 1 — 

£o , 


Due to the stopping criterion of Algorithm 1 we have for all i that for the last j —th iteration 
4 < Erprtj- Therefore, for all i large enough, due to (76) we get 


IPu - Pv ex 

\\JT U. -IT V 2 || C2 


(A) < 


and 


which implies that 


Pu- (id-P)vV 


0 -i) _ - 0 ,at)i 


A (A) - e 


9ej 2 0e 7 

< -4 + —4 


(Puk -(id-P)vV 


0-i) -0 ,at) 


— v,' 


< 0e,-. 


(79) 

(80) 


lb 2 (A) - 

Due to (Al), (80), and || • ||^» being a quasi-norm, it follows that v4 : = COARSE[(l — 0)ej,v^’ K ^ 
for i large enough, satisfies 


II^IU* < C 2 ( S )||(Puk.-(id-P)vp-^| 


< 


Ci(s)C 2 (s) (||(Pu)eJU* + || id 




so by (77) and ej = 3p K /Oej-i (see Algorithm 1), 


A'-'w^'Wa- 


< 


em# + Ci(»)ft(»)ii id-p|i(3^/e) J /— 1 ( e ->r'i 




We can conclude that for K > 0 large enough, and by the assumption (75), the solutions of the 
homogeneous part of this recursion converge to zero, and so 


e f- 1 iivp ) i u- < mil:, 


(81) 


uniformly with respect to j. And, due to (Al) we also have 

#supp(vp } ) < T(#supp(v ? 0) )) 

(Pu) e ,-(id-P)*i | U 


< £ -i/*n ro-n\ ha r>'w;W“i) n i/ 5 


< 


s/s 


(Pfl)«ilU* + l|id-P||||vp" 1) |Ui 


l/s 
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Therefore, by (77) and (81) we get with s/s = 1 + 5 that 


#supp(vp°) < T(#supp(v^)) < 


(i)> 


s/s 


s/s 


U 


(82) 


Next, recall that 0 < s < s and by Algorithm 1 we have £,; < ej for all j. Then (81) implies that for 
all j and i large enough we have e^/ s 1 ||v^ ||^s < ||u||^/ 


■ s / s 1m ”^" -- < ||u||Y/. In particular, for i large enough, 


4 /s 1 l|v»iu» < 


s/s 


(83) 


By the same argument, (82) yields for i large enough 


#supp(vi) < T(# supp(vi)) < el 


u 


(84) 


Note that the stopping criterion in Algorithm 2 implies that for large enough i we have e < e*. 
Therefore, from (84) we also get that 


1-1-<5 1+«S 

#supp(u e ) < T(#supp(u e )) <£ “||u|| A . 


(85) 


To prove part c) it is sufficient to show that the number of algebraic operations needed for each 
iterations, i.e., for the computation of v.j := SOLVER, A, f — A 1 (v,;_ 1 )], is < e “l 1+ ' 5)/s lleJlK 1+ ' 5 l/ s 


u 


A 3 


for all i large enough. Note that for small i the number of algebraic operations needed is bounded 
by C£ -( 1 + 5 )/s ||u||^ +<5)/s for £ —> 0. 

Due to L < 1 and by (67), there exists N £ N such that L N r* < e. Therefore, the N— th, 

i°g e -i+*)" 

-iog^oo- 15i ^ rw 


exit iteration, satisfies N := 


log -i (L) 

°. 

estimate the total number of operations by 


log i(£/r*) 


< 


. This implies that we can 


log _i(r*) 

-iog ‘(6)- 


< 


E 

i =0 


_ —(l+<5 )/s ||--^|i(l +^)/ s 


| (!+$)/ s 


log ^_ 1 (r*) 

-log - 1 (e)- lo /°_ l(L) 
0 £ o 


E 

i=0 




I (1+5)/s 


1 — f £, 


(+ I+i,/ ‘) 


log _i(r*) 

[-log i(e)~ lo /°_ l(L ) ] + l 


1 — £, 


-(l+5)/ S 

0 


< -—(1+5)/s I 


I (1+5)/* 


The last inequality is due to the fact that here we consider the asymptotic behavior for £ —» 0. 

Therefore to conclude the proof, it is sufficient to estimate the complexity of SOLVE. To do 
so one can follow the argument used in the proof of [40, Theorem 3.12 (II)] replacing [40, (3.19)] 
by (81) and using the assumptions (A1)-(A3) when relevant. Due to assumption (A3) one has 
r(#supp(g l °_ ) 1 )) < e“ 1/s ||v i _ 1 ||^ 4 / / and ||g,-i ) 1 ||^- < (1 + ||v»_iH^a). Therefore by assumption (A2) 
one has T(#supp(^\)) < e“ 1/;s ||v i _i||^ / / and ||fElU 5 ^ (l + ||vi-i|U»). By induction assumption 

we have that ||vp _1) ||^s < (l+||v,_i||_ 4 s) and T(# supp(vp _1) )) < e~ 1/s ||v^_i ||^/, that are trivially 
valid for j = 1. Therefore, again by (A2) one has 

HvPlU^a + llv.-ilU,), (86) 

and 

T(#supp(vp’ fe) )) < e~ 1/s ||v J _i||^ 4 / /, (87) 
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for all 0 < k < K. Therefore, by (A2) and (A3), the number of algebraic operations required to 
compute vp’ A ^ is < e~ 1 ^ s (l+ ||vj_ 1 ||^«) 1 / s + T(#supp(vj_ 1 )). Finally by assumption (Al) the ap¬ 
plication of v[ j) := COARSE[(l - 0) £j , v\ j ’ K) ] costs < T(# supp(vp’ A) )) + o(e" 1/s ~||vf ,K) ||^f) < 
e~ 1 ^ s (l + Hvi-ill^*) 1 /®. Thus, by (Al), (86) and (87) the induction hypothesis holds, i.e. ||v^|| A » < 
(1 + ||vi_i||_ 4 f) and T(#supp(vp ) )) < e“ 1/s ||v J _i||^ / / for any j G N 0 . This implies by induc¬ 
tion that computing from vp' -1 ^ takes a number of operations < + ||vj_iH^.) 1 /® + 

T(#supp(vj_i)). Since ej decreases geometrically and by formulas (83) and (84), one can estimate, 
as done before, the cost of the computation of v,; := SOLVER, A, 1 — Ai(v,;_i)] as a multiple of 

e i " 1/s ( 1 + l|v*-i|U»-) 1/s + 1 °g(£i)T(#supp(v i _i)) 

< e i /s (l + ||vi-i|U») 1/s + log(e i )e i _ 1 s ||u|| + 

S/a f f , \ 1 / s _l±l 1 + i 

£ i 3 [e? (i +||vi-i|U.-)J +\og(e i )e i _{ ||u|| + 

S/S S/s 1+5 1+5 

;$ £ i s l|u||+ + log(£i)e i _ 1 a ||u|| + 

1+5 1+5 

^ £i ‘ ll fl ll+ • 

In the last inequality we could ignore the log factor because of the arbitrary choice of 5 > 0, small 
as one wants. This concludes the proof. M 

REMARKS: 1. Theorem 6.1 ensures that our scheme is arbitrarily close to optimality, i.e., it is 
quasi-optimal; 

2. On the one hand, it is a very difficult theoretical problem to prove that P is bounded on A* for 
all t G (0, s) for wavelet frames. On the other hand, this condition has been verified numerically in 
[48, 14] in case of wavelet frame discretization, by observing the optimal convergence of SOLVE. 
According to [40, Remark 3.13], the boundedness of P on A t for all t G (0, s) is (almost) a necessary 
requirement for the scheme to behave optimally. 

There exists frames, for example time-frequency localized Gabor frames (and more generally all 
intrisically polynomially localized frames [25, 13]), for which the boundedness of the corresponding 
P has been proven rigorously, see [13, Theorem 7.1 in Section 7]. Therefore, for certain operator 
equations (for example, certain pseudodifferential equations appearing, e.g., in wireless communi¬ 
cation modeling [29]) the optimal application of SOLVE based on Gabor frame discretizations is 
justified theoretically. 

One can avoid requiring the boundedness of P and obtain a fully optimal scheme replacing 
SOLVE by its modified version modSOLVE (see [40]). One assumes then, without loss of gen¬ 
erality, that the number of algebraic operations needed to compute in (A3) is bounded by 
< Vs llv,]]^. The procedure modSOLVE requires the definition of an implementable alterna¬ 
tive projector P onto ran(A) and it is possible if a suitable wavelet frame expansion is constructed. 
In Section 7, we present such a construction that allows for a straightworward definition of the 
admissible P following the recipe in [40]. 

If IF is a R.iesz basis, then the scheme is certainly optimal and our result confirms the one in [10, 
Theorem 7.5]. 

3. Due to Theorem 5.1, (76) holds if the physical parameters of the MHD problem allow for 
choosing L, e 0 and r* such that 

0 < L =: r* 7 < £o < r* < y _1 

for 7 = 7 (.A, A,A{) \= 2||A _ 1 | ran ( A )||||A 1 ||||F|| 3 . In particular, if 0 < 7 < 1 is small enough then 
(76) is satisfied. We show next that, for the particular MHD equations and in case T is a Riesz basis, 
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the dependence of 7 on the viscosity 77 and the electric resistivity a 1 can be expressed explicitly. 
Note that 


<Au, u) = {{D^)~ 1 FAF*Dv 1 u,u) 

= (AF*D~ 1 u,F*D v 1 u) 

> a||F*£)^- 1 u||v 

= a il rXD?*) n fn\W ~ a ll-Py 1 u||vd 


^ 11 U 11 (At) • 


Therefore, if T is a Riesz basis, then ran(A) = I 2 (A/") and 

l|A|“n(A)ll = l!A- 1 ||<a- 1 , (88) 

where from (11) we have a := ^ ao — 2||ai|| ||(u 0 , Jo)||x (u and ao = c(fl) min{r/, cr —1 }. This 
implies that if 77 , <r ~ 1 are large enough then 7 < 1 can be made sufficiently small. The norm 
equivalence || ^neAx(-^ ) v 1 u)n/n|| v ~ H-D^ullvd va lid only if F is a Riesz basis and it does not 
hold for pure frames. It is possible to get an estimate similar to ( 88 ) under the additional assumption 
on the frame F that D^ 1 ran(A) C ran(F), which is to be verified for any particular frame under 
consideration. 


7 Construction of aggregated wavelet frames 

In this section we present the construction of suitable multiscale frames on bounded domains for 
the MHD problem described in Sections 2-3. They are in fact Gelfand frames for (V,7Y, V'), which 
ensures the correct application of Algorithms 1 and 2. In particular, we show that for the special 
case T = Ti = Si we can use the wavelet bases constructed using [46, Chapter 2, Definition 5] from 
a biorthogonal system on Z^fl) satisfying [46, Chapter 2, Assumption 3], where fl is some open and 
bounded domain, see for example [16, 17, 18]. 

We start by studying the properties of the functions in V. Using 

V := {(v, K) G X (ViJf) : b((v,K), (q, ip)) = 0 for all (qrf) G 

and setting 7 /) = 0 in the definition of the form b we obtain 

I (V • v)q = 0 for all q G ^(U). 

Jn 

By [7, Lemma l.a)] we get that V • v = 0. Thus, v G V(div; fi) satisfying v|r = 0, where 

V(div; fi) :={vGl 2 (!l): V • v = 0}. 

Therefore, v is in V(div;fl) fl Hq(U). A detailed discussion of construction of a wavelet basis for 
such a space is given, for example, in [44, 45, 46]. 

Note that the application of V to ip yields the same results regardless if ip G H 1 (fl) or if) G 
i7 1 (fl) \ R. Thus, substituting q = 0 into the definition of the form b we obtain 


/ K(Vi(>) = 0 for all if) G H\Q) \ R. (89) 

Jn 

Corollary 3.4. in [26] yields, for simply-connected U, that K = V x £ with £ G H(curl; O), where 
H(curl; fl) := {£ G L 2 , Vx(G ^ 2 }- There is also another way to describe the space, to which K 
belongs. 
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Proposition 7.1. Let K and ip be in £ 2 (ft) and F^ft), respectively. Then 

[ K(Vil’)= 0 for all ip G F x (ft) 

Jn 

if and only ifV-K = 0 and K • n|r = 0. 

Proof. ” => “ Integration by parts yields 

[ (V • K)ip = - [ K{Vip) + / K -nip for all ip G if 1 (ft) \ R. (90) 

do do dr 

And, by the assumption together with (89) we get that (90) reduces to 

f (V • K)ip = 0 for all ip G F^ft) C F^ft) \ R. 
do 

This implies that V • FT = 0 due to F _ 1 (ft) being a dual for Fg(ft). Thus, (90) becomes 

J K ■ nip = 0 for all ip G F^ft) \ R, 

which implies IT • n|r = 0 as IT • n|r is in F 1 / 2 (T)', the dual of F 1 / 2 (T), and ip G F 1 / 2 (T). The 
claim follows. ”<*=“ follows directly from (90). I 


In the notation of [46, (2.26)], by Proposition 7.1 K G Vo(div;ft) := Ho(div;ft) n V(div;ft), 
where by [26, Theorem 2.6] 

H 0 (div; ft) = {K G L 2 {Q) : V • K G £ 2 (ft), IT ■ n| r = 0}. 

Thus, for example, in case ft = (0, l) n , the existence of a Riesz basis for V 0 (div; ft) follows from [46, 
Theorem 10 and Proposition 5, p. 100]. Therefore, we are guaranteed that a wavelet frame needed 
for discretizing Problem 3 exists, at least, for any domain which is an affine image of (0,1)" [46, p. 
103]. Restricting our problem to 2D, we could even work with divergence-free wavelets on domains 
ft C R 2 which are conformal images of (0, l ) 2 [46, p. 104], 

We discuss next the methods for obtaining pure frames for V on more general domains. The 
following two approaches are of interest: 

• Modify the construction of pure frames for H s (ft) given in [13] and done using the ODD 

(Overlapping Domain Decomposition) technique. Using ODD assume that {ftj}^£ 1? M G N, 
are overlapping subdomains such that ft = u(]£ 1 fti. Such domains could be assumed to be affine 
images of the reference domain □ := (0, l) n (see Figure 4). Moreover, assume that for each 1 < 
i < M a divergence-free wavelet basis d', := is given for V(div; ft;) nHj(ftj). 

Show that d' := U^d/; = {ip)^}j>— i,kejy ,i=\,...,M is a Gelfand frame for V(div; ft) D H^ft). 

• Given a wavelet Gelfand frame for F 1 (ft), possibly constructed by ODD, follow a similar 
strategy as in [44, 45, 46] and construct a divergence-free Gelfand frame for V(div; ft)nHj(ft). 

These systems, called aggregated divergence-free wavelet frames, allow us to avoid dealing with 
interfacing patches used in disjoint domain decompositions (DDD) (see [18, 46]). DDD are usually 
rather complicated to implement and can yield ill-conditioned systems (see [46, p. 104, sec. More 
General Domains]). 

We continue by presenting some preliminary results on constructions of the desired pure wavelet 
frames described in the first approach above. Assume that C := {ft;}^ is an overlapping, relatively 
compact covering of ft such that 
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Figure 3: Example of an Overlapping Domain Decomposition of a polygonal domain in 2D by means 
of patches which are affine images of (0, l) n . 


(Cl) there exist affine (or even conformal if n = 2) maps /«$:□—> fij, for all? = 1,, M. 

The set of admissible domains fl is restricted by raising condition (Cl), e.g., the boundary of fl 
has to be piecewise smooth. Nevertheless, the particularly attractive case of polyhedral domains is 
still covered. Thus, the above assumptions on the parametrizations k, are by no means principal 
restrictions on the shape of fl. 

We assume that the wavelets we use below have sufficient regularity and number of vanishing mo¬ 
ments. We consider a template vector-valued wavelet divergence-free basis \F n = kej° 

for Vg := V(div;D) fl Hg(D), s > 0, and its biorthogonal dual \F n = i^k}j>j 0 keJ°' ^ 
that \F n is a frame for V^. Our aim is to show that the system 


vF : (V’j 1 fc)(j,j,fc)eA) 


where 




1 (x)) V)) 

| det VKj(K l “ 1 (a;)) | 1/2 


for all i = 1,..., M, j > jo, k G Jj x G f2*, 


(91) 


(92) 


and 4>j tk (x) = 0 for x € fl \ fli, is a frame for V s := V(div; fl) fl Hg(fl). Anagously, we define its 
local duals by 


V’UO) := 


^a:)) ^a:)) 


detVft^Kj x (x)) 


11/2 


for all i = 1,..., M, j > jo, k G J~, x G f \ 


(93) 


and ipj k (x) = 0 for a: G fl \ flj. 

First of all observe that under assumption (Cl) each fF* := (' , Pj tk )j> JU kej * 3 a S a i n a divergence- 

free wavelet frame for Vf := V(div;fli) nHg(fl;). Moreover, by [44, 45, 46] for Ui = JT fc c 7 ',fcV’},/ s 
the following norm equivalence holds 


1/2 


u * H» 


' 2 2sj I 




s > 0. 


(94) 


We next prove the following two auxiliary lemmas. 
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Lemma 7.2. Let tH 0,1 and \I / n ' 2 be two generic wavelet bases on □ with sufficient regularity and 
number of vanishing moments. We do not necessarily require that these wavelets vanish on the 
boundary. Moreover, let LI be affine (or conformal if n = 2) image of D. Let the corresponding 
wavelet systems IH 1 = (ifj k)j>j 0 keJ a<1 anc ^ ^ = ^jk)j>j 0 k&j D ’ 2 on ^ con structed by lifting 
as in (92). Then the cross Gramian matrix 

M* 1 ’* 2 := ((V’j.fe) V’|',fe'))(i,fe)eA 1 ,0' , ,fe')eA 2 ( 95 ) 

defines a bounded operator from i 22 <o(A 2 ) onto i 22 ^(A 1 ) with 




f \ 

1/2 

(-2,2 s i (A) : — < 

( C j,k)(j,k)£A '■ 

E 22sj \ c ^\ 2 

\j>jo,kejP ^ 

< 00 

> 


for any |s| < so, where s o G M+ depends on the regularity and the number of vanishing moments of 
and \F n ’ 2 . 

Proof. The result is well known and it can be immediately achieved, e.g., by combining [13, Theorem 
6.5] and [13, Proposition 8.2], where in the latter one replaces p\ with p 2 . S 


Lemma 7.3. Let C := {Lli}ff 1 be an overlapping, relatively compact covering of LI satisfying (Cl). 
Let the projections of L 2 (Ll) onto Vj := V(div, fi$) be 


p Vii u) := E ( u ’4>j, k Wj,k’ 

j>j 0 ,kejP 


(96) 


Then 


(i) 


( 


1/2 


< 


u I|h« for any i = 1 ,... ,M and all u G H s (fl); 


E 2 2 «|( u ,^ ifc )| 2 

\j>jo,k€jp ) 

(ii) ||TVi(u)|| h s < || u I|h* for any i = 1,... ,M and all u G H s (fl). 

Proof. Fix i = 1 Cleary u|o i G H s (f2j), if u G H s (0). It is well-known that u|n i = 

Lkij>j 0 ,kejf 


^2j>j g ke jO,« (u, i’j’Di’j k f° r some suitable wavelet basis {V’, k^x a -° constructed by lifting a 
given reference basis on □. Moreover, 


1/2 


E 2 2 ^|< uE ;°)| 2 


< 


u H« 


(97) 


\j>jo,k&J 




and 


j>jo,keJ p ° 

Combining (97) and Lemma 7.2 we get (i). Due to (94) we have 



f \ 

1/2 

f \ 

||7V i (u)|| H = ~ 

E 2 2 ^'|(u,4 ifc ,)| 2 

\i’>jo,k'ej° j 

< 

E 2 2 ^l(uE;°}| 2 

\j>30,kej°’° ) 


1/2 


< 


U H«. 
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We next present the frame construction that allows us to work with Q, which are not necessarily 
affine images of (0, l) n , and provides the characterization of V s for any s > 1. We require that C is 
such that 

(C2) it holds 

^V«( v )ln i \uj“ =2 n ii . = v lni\u* 2 n ifc > 

for any v G V s such that v| asupp(v) = 0, C supp(v) C for {n,C {1, 

without loss of generality i :=i\. 

REMARKS: 

1. Note that the domain decomposition satisfying (C2) is necessarily overlapping. If it were not 

the case, then it would hold that f \ = fl,; \ and, by (98), Py 4 v = v|n 4 . This, would 

imply that v|an« = Pyvlac^ = 0, which is not necessarily the case for an arbitrary v G V s . 

2. To illustrate our frame construction we consider Q = fli U CI 2 , which is not an affine image 
of (0, l) 2 , as in Figure 5. It holds Py^v = v on Q \ LI 2 and Pyv G V s , s > 1. This, implies that 
v — Piyv G V s . Note that by construction supp(v — P^v) C fi 2 and v — Py 1 v|on 2 = 0- Therefore, 
we obtain the following decomposition v = Py x v + Py 2 (v — Pyv). 


/ 



/ 





/ 

/ Q.2 

/ 

/ 

/ 


Figure 4: An example of Overlapping Domain Decomposition satisfying (Cl) and (C2). 

The following Theorem generalizes the argument in Remark 2. above to domain decompositions 
consisting of more than two patches. 

Theorem 7.4. Let C := be an overlapping, relatively compact covering of D satisfying 

(C1)-(C2). Let T := (V’) fe)(i,i,fe)eA be as in (91) and Pv t as in Lemma 7.3. For u G V s , s > 1, 
consider the functions and domains defined by the following recursion 


u<°> = u, 


= n. 


Then 


Un+i = Py„ +1 u(") 

U (n+1) = u (n) _ Un+1 5 n = 0,..., M - 1. 
f 2 ( n+1 ) = supp(u(")) 


(99) 


M 


M 


u = E Ui = E El < u(i 


2=1 






□ 


( 100 ) 


with 


U H a ~ 


M 


( 

E E 

\* =1 3>3o,kejp 


1/2 




( 101 ) 


In particular, 4/ is an L 2 frame for Vo(div; D) and a Gelfand frame for (V s , Vo(div; fi), (V s )'). 
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Proof. By property (C2), Lemma 7.3 and by induction we get that u*") € V s with ||u(")||h s < 
||u||h s , u^"^|af 2 (' i + 1 ) = 0 and C U f£ n+1 fli. Thus, fl( M ) = 0 and, therefore, in the recursive 
definition (99) we have 0 < n < M — 1. This implies also that supp(u^ M_1 l) C VLm and u m = 
Pv M vf M ~ 1 ' > = And by induction we gets = u —Uj. This immediately implies 


M—l / M—l \ M-1 M 

u = + u ~ yy u, j = yy Ui +u ( m_i) = yy u 4 

i= 1 V i=1 / i=l i=l 

and, consequently, the wavelet decomposition of u 


M 


M 


M 


" = £”< = E / v."““ I> = £ E 


i=i i=i *=i 

By an application of Lemma 7.3 (i) and property (94) one has 

/ 


1/2 


u||h* ^ 


M 


£ yy 

\ i=1 j>jo,kejp 


Conversely, by Lemma 7.3 (i) 

( 


1/2 


M 


M 


1/2 


yy yy 2 2 «i<u< < - i ),^*. jfc >i 2 1 < (y^ nu^- 1 ) 

\ i=1 i>jo,kejp 


— 1) l|2 

H s 


< llul 


H s • 


Note that the norm equivalences above hold for also for s = 0 and, thus, by density of V s in 
V 0 (div; fi), \I> is an L 2 frame for V 0 (div; 0) and a Gelfand frame for (V s , V 0 (div; O), (V s )'). I 


REMARKS: 1. As mentioned in Remark 2. after Theorem 6.1, the wavelet frame expansion (100) 
allows for the construction of an admissible projector P for the optimal application of the procedure 
modSOLVE. Such construction follows the one in [40, Section 4.4]. 

2. The aggregate wavelet frames constructed in Theorem 7.4 are multiscale and produce co¬ 
efficients that can be naturally structured into suitable trees, as the ones in [17, 18]. Thus, the 
algorithms defined in [11], see (A3), for the evaluation of nonlinear functionals can also be applied in 
our case. The optimality also requires the characterization of Besov spaces in terms of norm equiv¬ 
alences. As soon as the local divergence-free bases on □ provides such a characterization, one can 
show similarly to Theorem 7.4 that the global divergence-free basis provides the characterization of 
Besov spaces on the whole domain. We postpone the details to a forthcoming paper. 

3. In practice, one never uses the full frame ’S' := but rather 

i&j := {'(/>} k}-i<j<J,keji,i=i,...,M (f° r example, in the concrete implementation of the schemes dis¬ 
cussed in [48], an upper bound in terms of concretely used scale levels j must be enforced for compu¬ 
tational (memory/storage) limits, without spoiling the accuracy). Therefore, the first approach men¬ 
tioned above is “ready-to-use”. It is well-known that any finite system is already a (Gelfand) frame 
for its span. Using this fact one can consider solving Problem 4 in span(4/ j) C V(div; f l) n Hq(O) 
for J > 0 large. The scheme we propose is well-defined in this case. 


8 Conclusion 

We present the first convergent and implementable adaptive scheme based on frame decompositions 
for the numerical integration of magnetohydrodynamic flows. Certainly, it is impossible to construct 
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a universal scheme suitable for any physical problem. Thus, the convergence of the algorithm is 
ensured under certain assumptions on the physical parameters. Such assumptions are standard 
when dealing with nonlinear problems of fluid dynamics and magnetohydrodynamics and still lead 
to the analysis of physically realistic problems, see [34, 36, 35]. 

We show that under Dirichlet boundary conditions on the velocity and the zero flux boundary 
conditions on the electric current the velocity and electric current are divergence-free vector valued 
functions. The discretization of the problem is, therefore, realized by expanding the solution with 
respect to certain divergence-free wavelet frames constructed on Overlapping Domain Decomposi¬ 
tions. The construction based on ODD is definitively a breakthrough in comparison with [44, 45, 46]. 
There divergence-free wavelet bases have been defined in a relatively simple way essentially only 
for domains that are affine images of cubes. Our construction of frames allows for rather general 
polygonal domains and avoids the use of fictitious domain techniques, see, e.g., [15]. These frames 
can be also used for the discretization of Stokes and Navier-Stokes equations [19]. Of importance 
is that these aggregate divergence-free wavelet frames preserve multiscale properties of standard 
wavelet bases and their capability to characterize certain function spaces. These properties ensure 
the optimal convergence rates and complexity of the scheme we propose, if the solutions are in 
certain sparseness classes of functions. It has not yet been theoretically proven that the solutions 
belong to such function spaces (e.g., Besov spaces), numerical simulations [1] though support this 
assumption. 
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